/*
 * IlutpPreconditioner.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
using System.Collections;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// This class performs an Incomplete LU factorization with drop tolerance
    /// and partial pivoting. The drop tolerance indicates which additional entries
    /// will be dropped from the factorized LU matrices.
    /// </summary>
    /// <remarks>
    /// The ILUTP-Mem algorithm was taken from: <br/>
    /// ILUTP_Mem: a Space-Efficient Incomplete LU Preconditioner
    /// <br/>
    /// Tzu-Yi Chen, Department of Mathematics and Computer Science, <br/>
    /// Pomona College, Claremont CA 91711, USA <br/>
    /// Published in: <br/>
    /// Lecture Notes in Computer Science <br/>
    /// Volume 3046 / 2004 <br/>
    /// pp. 20 - 28 <br/>
    /// Algorithm is described in Section 2, page 22
    /// </remarks>
    /// <include file='../../../../../examples/trunk/LinearAlgebra/Solvers/Preconditioners/IlutpPreconditioner.xml' path='example'/> 
    public sealed class IlutpPreconditioner : IPreconditioner
    {
        /// <summary>
        /// The original coefficient matrix.
        /// </summary>
        private SparseMatrix mCoefficientMatrix;
        /// <summary>
        /// The decomposed upper triangular matrix.
        /// </summary>
        private SparseMatrix mU;
        /// <summary>
        /// The decomposed lower triangular matrix.
        /// </summary>
        private SparseMatrix mL;

        /// <summary>
        /// A flag used to indicate if the preconditioner has been created, i.e.
        /// if the <b>mLuDecomposition</b> matrix has been filled with the correct
        /// values.
        /// </summary>
        private bool mWasPreconditionerCreated;

        /// <summary>
        /// The array containing the pivot values.
        /// </summary>
        private int[] mPivots;
        /// <summary>
        /// The fill level.
        /// </summary>
        private double mFillLevel = 200.0;
        /// <summary>
        /// The drop tolerance.
        /// </summary>
        private double mDropTolerance = 0.0001;
        /// <summary>
        /// The pivot tolerance.
        /// </summary>
        private double mPivotTolerance;

        /// <summary>
        /// Creates a new <b>IlutpPreconditioner</b> instance.
        /// </summary>
        /// <param name="coefficientMatrix">
        /// The <see cref="Matrix"/> upon which this preconditioner is based. Note that the 
        /// constructor takes a general matrix type. However internally the data is stored 
        /// as a sparse matrix. Therefore it is not recommended to pass a dense matrix to the
        /// constructor.
        /// </param>
        /// <exception cref="ArgumentNullException">If <paramref name="coefficientMatrix"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="coefficientMatrix"/> is not a square matrix.</exception>
        public IlutpPreconditioner(Matrix coefficientMatrix)
        {
            if (coefficientMatrix == null)
            {
                throw new ArgumentNullException("coefficientMatrix", Strings.NullParameterException);
            }

            if (coefficientMatrix.Rows != coefficientMatrix.Columns)
            {
                throw new NotSquareMatrixException();
            }

            mCoefficientMatrix = new SparseMatrix(coefficientMatrix);
        }

        /// <summary>
        /// Returns the decomposed upper triangular matrix.
        /// </summary>
        /// <remarks>
        /// This method is used for debugging purposes only and 
        /// should normally not be used.
        /// </remarks>
        /// <returns>The upper triangular matrix.</returns>
        internal SparseMatrix UpperTriagonalMatrix()
        {
            return mU;
        }

        /// <summary>
        /// Returns the decomposed lower triangular matrix.
        /// </summary>
        /// <remarks>
        /// This method is used for debugging purposes only and 
        /// should normally not be used.
        /// </remarks>
        /// <returns>The lower triangular matrix.</returns>
        internal SparseMatrix LowerTriagonalMatrix()
        {
            return mL;
        }

        /// <summary>
        /// Returns the pivot array. This array is not needed for normal use because
        /// the preconditioner will return the solution vector values in the proper order.
        /// </summary>
        /// <remarks>
        /// This method is used for debugging purposes only and 
        /// should normally not be used.
        /// </remarks>
        /// <returns>The pivot array.</returns>
        internal int[] Pivots()
        {
            int[] result = new int[mPivots.Length];
            Buffer.BlockCopy(mPivots, 0, result, 0, mPivots.Length * Constants.SizeOfInt);
            return result;
        }

        /// <summary>
        /// Gets or sets the amount of fill that is allowed in the matrix. The
        /// value is a fraction of the number of non-zero entries in the original
        /// matrix. The standard value is 200.
        /// </summary>
        /// <remarks>
        /// <para>
        /// Values should always be positive and can be higher than 1.0. A value lower
        /// than 1.0 means that the eventual preconditioner matrix will have fewer
        /// non-zero entries as the original matrix. A value higher than 1.0 means that
        /// the eventual preconditioner can have more non-zero values than the original 
        /// matrix.
        /// </para>
        /// <para>
        /// Note that any changes to the <b>FillLevel</b> after creating the preconditioner
        /// will invalidate the created preconditioner and will require a re-initialization of
        /// the preconditioner.
        /// </para>
        /// </remarks>
        public double FillLevel
        {
            get
            {
                return mFillLevel;
            }
            set
            {
                mFillLevel = value;
                // We have changes here so we'll need to recreate the
                // preconditioner.
                mWasPreconditionerCreated = false;
            }
        }

        /// <summary>
        /// Gets or sets the absolute drop tolerance which indicates below what absolute value
        /// an entry will be dropped from the matrix. The standard value is 0.0001.
        /// </summary>
        /// <remarks>
        /// <para>
        /// The values should always be positive and can be larger than 1.0. A low value will
        /// keep more small numbers in the preconditioner matrix. A high value will remove 
        /// more small numbers from the preconditioner matrix.
        /// </para>
        /// <para>
        /// Note that any changes to the <b>DropTolerance</b> after creating the preconditioner
        /// will invalidate the created preconditioner and will require a re-initialization of
        /// the preconditioner.
        /// </para>
        /// </remarks>
        public double DropTolerance
        {
            get
            {
                return mDropTolerance;
            }
            set
            {
                mDropTolerance = value;
                // We have changes here so we'll need to recreate the
                // preconditioner.
                mWasPreconditionerCreated = false;
            }
        }

        /// <summary>
        /// Gets or sets the pivot tolerance which indicates at what level pivoting will
        /// take place. The standard value is 0.0 which means pivoting will never take place.
        /// </summary>
        /// <remarks>
        /// <para>
        /// The pivot tolerance is used to calculate if pivoting is necessary. Pivoting
        /// will take place if any of the values in a row is bigger than the 
        /// diagonal value of that row divided by the the pivot tolerance, i.e. pivoting
        /// will take place if <b>row(i,j) > row(i,i) / PivotTolerance</b> for
        /// any <b>j</b> that is not equal to <b>i</b>.
        /// </para>
        /// <para>
        /// Note that any changes to the <b>PivotTolerance</b> after creating the preconditioner
        /// will invalidate the created preconditioner and will require a re-initialization of
        /// the preconditioner.
        /// </para>
        /// </remarks>
        public double PivotTolerance
        {
            get
            {
                return mPivotTolerance;
            }
            set
            {
                mPivotTolerance = value;
                // We have changes here so we'll need to recreate the
                // preconditioner.
                mWasPreconditionerCreated = false;
            }
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        /// <remarks>
        /// <para>
        /// When called this method always performs the creation of the preconditioner, even
        /// if the preconditioner has been created previously and no changes have taken place.
        /// </para>
        /// <para>
        /// Note that the created preconditioner will be invalid if changes to the matrix 
        /// occur after calling <b>CreatePreconditioner</b>.
        /// </para>
        /// </remarks>
        public void CreatePreconditioner()
        {
            // The creation of the preconditioner follows the following algorithm.
            //spaceLeft = lfilNnz * nnz(A)
            //for i = 1, .. , n
            //{
            //    w = a(i,*)
            //    for j = 1, .. , i - 1
            //    {
            //        if (w(j) != 0)
            //        {
            //            w(j) = w(j) / a(j,j)
            //            if (w(j) < dropTol)
            //            {
            //                w(j) = 0;
            //            }
            //            if (w(j) != 0)
            //            {
            //                w = w - w(j) * U(j,*)
            //            }
            //        }
            //    }
            //
            //    for j = i, .. ,n
            //    {
            //        if w(j) <= dropTol * ||A(i,*)||
            //        {
            //            w(j) = 0
            //        }
            //    }
            //
            //    spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
            //    lfil = spaceRow / 2  // space for this row of L
            //    l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
            //
            //    lfil = spaceRow - nnz(L(i,:))  // space for this row of U
            //    u(i,j) = w(j) for j = i, .. , n // only the largest lfil - 1 elements
            //    w = 0
            //    
            //    if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
            //    {
            //        pivot by swapping the max and the diagonal entries
            //        Update L, U
            //        Update P
            //    }
            //    spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
            //}

            mL = new SparseMatrix(mCoefficientMatrix.Rows);
            // Create the upper triangular matrix and copy the values
            mU = new SparseMatrix(mCoefficientMatrix.Rows);
            // Create the pivot array
            mPivots = new int[mCoefficientMatrix.Rows];
            for (int i = 0; i < mPivots.Length; i++)
            {
                mPivots[i] = i;
            }
            Vector workVector = VectorBuilder.CreateVector(mCoefficientMatrix.Rows, VectorType.Dense);
            Vector rowVector = VectorBuilder.CreateVector(mCoefficientMatrix.Columns, VectorType.Dense);
            int[] indexSorting = new int[mCoefficientMatrix.Rows];
            // spaceLeft = lfilNnz * nnz(A)
            int spaceLeft = (int)mFillLevel * mCoefficientMatrix.NonZeros;
            // for i = 1, .. , n
            for (int i = 0; i < mCoefficientMatrix.Rows; i++)
            {
                // w = a(i,*)
                mCoefficientMatrix.Row(i, workVector);
                // pivot the row
                PivotRow(workVector);
                double vectorNorm = workVector.InfinityNorm();
                // for j = 1, .. , i - 1
                for (int j = 0; j < i; j++)
                {
                    // if (w(j) != 0)
                    // {
                    //     w(j) = w(j) / a(j,j)
                    //     if (w(j) < dropTol)
                    //     {
                    //         w(j) = 0;
                    //     }
                    //     if (w(j) != 0)
                    //     {
                    //         w = w - w(j) * U(j,*)
                    //     }
                    if (workVector[j] != 0.0)
                    {
                        // Calculate the multiplication factors that go into 
                        // the L matrix
                        workVector[j] = workVector[j] / mU[j, j];
                        if (System.Math.Abs(workVector[j]) < mDropTolerance)
                        {
                            workVector[j] = 0.0;
                        }
                        // Calculate the addition factor
                        if (workVector[j] != 0.0)
                        {
                            // vector update all in one go
                            mU.Row(j, rowVector);
                            // zero out columnVector[k] because we don't need that
                            // one anymore for k = 0 to k = j
                            for (int k = 0; k <= j; k++)
                            {
                                rowVector[k] = 0.0;
                            }
                            rowVector.Multiply(workVector[j]);
                            workVector.Subtract(rowVector);
                        }
                    }
                }
                //    for j = i, .. ,n
                for (int j = i; j < mCoefficientMatrix.Rows; j++)
                {
                    // if w(j) <= dropTol * ||A(i,*)||
                    // {
                    //     w(j) = 0
                    // }
                    if (System.Math.Abs(workVector[j]) <= mDropTolerance * vectorNorm)
                    {
                        workVector[j] = 0.0;
                    }
                }
                // spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
                int spaceRow = spaceLeft / (mCoefficientMatrix.Rows - i + 1);
                // lfil = spaceRow / 2  // space for this row of L
                int fillLevel = spaceRow / 2;
                FindLargestItems(0, i - 1, indexSorting, workVector);
                // l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
                int lowerNonZeroCount = 0;
                int count = 0;
                for (int j = 0; j < i; j++)
                {
                    if ((count > fillLevel) || (indexSorting[j] == -1))
                    {
                        break;
                    }
                    mL[i, indexSorting[j]] = workVector[indexSorting[j]];
                    count += 1;
                    lowerNonZeroCount += 1;
                }
                FindLargestItems(i + 1, mCoefficientMatrix.Rows - 1, indexSorting, workVector);
                // lfil = spaceRow - nnz(L(i,:))  // space for this row of U
                fillLevel = spaceRow - lowerNonZeroCount;
                // u(i,j) = w(j) for j = i + 1, .. , n // only the largest lfil - 1 elements
                int upperNonZeroCount = 0;
                count = 0;
                for (int j = 0; j < mCoefficientMatrix.Rows - i; j++)
                {
                    if ((count > fillLevel - 1) || (indexSorting[j] == -1))
                    {
                        break;
                    }
                    mU[i, indexSorting[j]] = workVector[indexSorting[j]];
                    count += 1;
                    upperNonZeroCount += 1;
                }
                // Simply copy the diagonal element. Next step is to 
                // see if we pivot
                mU[i, i] = workVector[i];
                // if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
                // {
                //     pivot by swapping the max and the diagonal entries
                //     Update L, U
                //     Update P
                // }
                // Check if we really need to pivot. If (i+1) >=(mCoefficientMatrix.Rows -1) then
                // we are working on the last row. That means that there is only one number
                // And pivoting is useless. Also the indexSorting array will only contain
                // -1 values.
                if ((i + 1) < (mCoefficientMatrix.Rows - 1))
                {
                    if (System.Math.Abs(workVector[i]) < mPivotTolerance * System.Math.Abs(workVector[indexSorting[0]]))
                    {
                        // swap columns of u (which holds the values of A in the
                        // sections that haven't been partitioned yet.
                        SwapColumns(mU, i, indexSorting[0]);
                        //  Update P
                        int temp = mPivots[i];
                        mPivots[i] = mPivots[indexSorting[0]];
                        mPivots[indexSorting[0]] = temp;
                    }
                }
                // spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
                spaceLeft -= (lowerNonZeroCount + upperNonZeroCount);
            }
            for (int i = 0; i < mL.Rows; i++)
            {
                mL[i, i] = 1.0;
            }
            mWasPreconditionerCreated = true;
        }

        private void PivotRow(Vector row)
        {
            Hashtable knownPivots = new Hashtable();
            // pivot the row
            for (int i = 0; i < row.Count; i++)
            {
                if ((mPivots[i] != i) && (!PivotMapFound(knownPivots, i)))
                {
                    // store the pivots in the hashtable
                    knownPivots.Add(mPivots[i], i);

                    double t = row[i];
                    row[i] = row[mPivots[i]];
                    row[mPivots[i]] = t;
                }
            }
        }

        private bool PivotMapFound(Hashtable knownPivots, int currentItem)
        {
            if (knownPivots.ContainsKey(mPivots[currentItem]))
            {
                if (knownPivots[mPivots[currentItem]].Equals(currentItem))
                {
                    return true;
                }
            }

            if ((knownPivots.ContainsKey(currentItem)))
            {
                if (knownPivots[currentItem].Equals(mPivots[currentItem]))
                {
                    return true;
                }
            }

            return false;
        }

        private static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
        {
            for (int i = 0; i < matrix.Rows; i++)
            {
                double temp = matrix[i, firstColumn];
                matrix[i, firstColumn] = matrix[i, secondColumn];
                matrix[i, secondColumn] = temp;
            }
        }

        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
        {
            // Copy the indices for the values into the array
            for (int i = 0; i < upperBound + 1 - lowerBound; i++)
            {
                sortedIndices[i] = lowerBound + i;
            }
            for (int i = upperBound + 1 - lowerBound; i < sortedIndices.Length; i++)
            {
                sortedIndices[i] = -1;
            }
            // Sort the first set of items.
            // Sorting starts at index 0 because the index array
            // starts at zero
            // and ends at index upperBound - lowerBound
            IlutpHeapSorter.SortDoubleIndicesDecreasing(0, upperBound - lowerBound, sortedIndices, values);
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <remarks>
        /// If the preconditioner was not created before calling the <b>Solve</b> method it will be
        /// called by the <b>Solve</b> method. The <b>CreatePreconditioner</b> is only called once
        /// even if the <b>Solve</b> method is called multiple times. Note that changes to the
        /// <b>PivotTolerance</b>, <b>DropTolerance</b> or the <b>FillLevel</b> will trigger a 
        /// rebuild of the preconditioner.
        /// </remarks>
        /// <param name="vector">The right hand side vector.</param>
        /// <param name="result">The left hand side vector.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If <paramref name="vector"/> and <paramref name="result"/> do not
        /// have the same size.</para>
        /// <para> - or -</para>
        /// <para>If the size of <paramref name="result"/> is different the number of rows of the coefficient matrix.</para></exception>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the size of <paramref name="vector"/> is different the number of rows of the coefficient matrix.</para></exception>
        public void Solve(Vector vector, Vector result)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Count != vector.Count) || (result.Count != mCoefficientMatrix.Rows))
            {
                throw new NotConformableException(Strings.ParametersNotConformable);
            }

            if (!mWasPreconditionerCreated)
            {
                CreatePreconditioner();
            }
            // Solve equation here
            // Pivot(vector, result);
            // Solve L*Y = B(piv,:)
            Vector rowValues = VectorBuilder.CreateVector(mL.Rows, VectorType.Dense);
            for (int i = 0; i < mL.Rows; i++)
            {
                mL.Row(i, rowValues);

                double sum = 0.0;
                for (int j = 0; j < i; j++)
                {
                    sum += rowValues[j] * result[j];
                }

                result[i] = (vector[i] - sum);
            }

            // Solve U*X = Y;
            for (int i = mU.Rows - 1; i > -1; i--)
            {
                mU.Row(i, rowValues);

                double sum = 0.0;
                for (int j = mU.Rows - 1; j > i; j--)
                {
                    sum += rowValues[j] * result[j];
                }

                result[i] = 1 / rowValues[i] * (result[i] - sum);
            }

            // We have a column pivot so we only need to pivot the
            // end result not the incomming right hand side vector
            Vector temp = result.Clone();
            Pivot(temp, result);
        }

        private void Pivot(Vector vector, Vector result)
        {
            for (int i = 0; i < mPivots.Length; i++)
            {
                result[i] = vector[mPivots[i]];
            }
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <remarks>
        /// When using the <b>Solve</b> method consider carefully if the
        /// left hand side vector can be reused. If so it would be better to 
        /// use the overloaded <b>Solve</b> method.
        /// </remarks>
        /// <param name="vector">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        public Vector Solve(Vector vector)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (vector.Count != mCoefficientMatrix.Columns)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            Solve(vector, result);
            return result;
        }
    }
}
